<?php
class CStat {

	const STEP = 500;

	// Gamma function coef (Lanczos method)
	public static $GP = array(
		0.99999999999980993, 
		676.5203681218851,
		-1259.1392167224028,
		771.32342877765313, 
		-176.61502916214059, 
		12.507343278686905,
		-0.13857109526572012, 
		9.9843695780195716e-6, 
		1.5056327351493116e-7
	);

  public function g($t, $x) {
    return pow($t, $x - 1.0) * exp(-$t);
  }

  public function gammaLow($a,$x) {

    $dx = $x / (double)self::STEP;
    $sx = (double)0.0;
    $s = (double)0.0;

    for ($i = 0; $i < self::STEP; $i++) {
      $s += ($this->g($sx, $a) + $this->g($sx + $dx, $a)) / 2.0;
	  echo $s . '<br>';
      $sx += $dx;
    }
    return $s * $dx;

  }

  public function gammaHigh($a, $x) {
    return $this->gamma($a) - $this->gammaLow($a, $x);
  }

  public function gamma($x) {

    // Reflection formula
    if ($x < 0.5) {
      return PI() / (sin(PI() * $x) * gamma(1.0 - $x));
    } else {
      $x -= 1.0;
      $y = self::$GP[0];
      for ($i = 1; $i < count(self::$GP); $i++) {
        $y += self::$GP[$i] / ($x + (double)$i);
      }
      $t = $x + (double)7.5;
      return sqrt(2 * PI()) * pow($t, $x + 0.5) * exp(-$t) * $y;
    }

  }
}
?>